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Eduardo Mendez-Villuendas,^ Ivan Saika-Voivod*,^ and Richard K. Bowles^^ 

Department of Chemistry, University of Saskatchewan, Saskatoon, Saskatchewan, SIN 5C9, Canada 

(Dated: February 1, 2008) 

We examine the metastable liquid phase of a supercooled gold nanocluster by studying the free 
energy landscape using the largest solid-like embryo as an order parameter. Just below freezing, 
the free energy exhibits a local minimum at small embryo sizes and a maximum at larger embryo 
sizes which denotes the critical embryo size. At T = 660K the free energy becomes a monotonically 
decreasing function of the order parameter as the liquid phase becomes unstable, indicating we have 
reached a spinodal. In contrast to the usual mean-field theory predictions, the size of the critical 
embryo remains finite as the spinodal is approached. We also calculate the rate of nucleation, 
independently from our free energy calculations, and observe a rapid increase in its temperature 
dependence when the free energy barrier is in the order of kT. This supports the idea that freezing 
becomes a barrierless process around the spinodal temperature. 



PACS numbers: 61.46.Df, 64.60.Qb, 64.60.My 

When a liquid is cooled below its freezing temperature 
we generally expect the system to crystallize. However, 
the nucleation process requires the formation of a small 
embryo of the new stable phase that introduces an ener- 
getically unfavourable liquid-solid interface and creates a 
free energy barrier between the liquid and solid phases. 
As long as the fluctuations in the liquid only result in 
the formation of embryos smaller than the critical size 
needed to overcome the barrier, the system will remain 
fluid. This is the metastable liquid [ij. 

As the liquid is cooled further, the free energy bar- 
rier decreases in height, making nucleation more likely 
and shortening the lifetime of the metastable liquid. The 
question then arises: Is there a temperature below which 
the metastable liquid becomes unstable with respect to 
all fluctuations? Mean-field theories, such as the gradi- 
ent theory developed by Cahn and Hilliard [2], predict 
such a limit of stability, or spinodal, for first order phase 
transitions like the condensation of a gas or liquid-liquid 
phase separation in a mixture. They also predict that 
the size of the critical nucleus diverges as the spinodal 
is approached as a result of the divergence in the mean- 
fleld isothermal compressibility of the fluid 3], and that 
the nucleation lifetime should diverge, despite the fact 
that the free energy barrier is in the order of kT , as 
the dynamics become increasingly cooperative Q- How- 
ever, recent experiments of phase separating polymers 
and simulations of the Ising model suggest that the size 
of the critical embryo remains finite as the spinodal is 
approached [5]. 

Whether a deeply supercooled single component liq- 
uid exhibits a spinodal singularity with respect to the 
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crystal remains an open question Trudu et al Q 
studied freezing in a bulk Lennard- Jones fluid and found 
nucleation becomes a spatially diffuse and collective phe- 
nomenon when the system is deeply supercooled and sug- 
gested this was indicative of the presence of a mean-field 
spinodal. Recent nucleation experiments on water show 
nucleation times become extremely short when the liq- 
uid is highly compressed, thus defining a practical limit 
of stability to the liquid state 0] . These results provide 
strong but indirect evidence for the existence of a ther- 
modynamic limit of stability for the supercooled liquid 
state. 

In this letter, we directly locate the limit of stabil- 
ity of the liquid phase by calculating the free energy of 
the cluster using the largest-sized solid embryo as an or- 
der parameter. At temperatures just below freezing, the 
free energy exhibits a local minimum associated with the 
metastable liquid and a free energy barrier that separates 
this liquid from the solid phase. The height of the barrier 
decreases as the temperature is lowered and eventually 
disappears so that the free energy becomes a monotoni- 
cally decreasing function of the order parameter and the 
liquid phase becomes unstable. This provides the first 
direct measurement of the spinodal in a simple liquid 
system. 

A rigorous molecular theory of a metastable system 
requires the introduction of constraints that prevent the 
system from accessing regions of phase space that will 
cause the system to evolve towards the more stable state. 
In the context of a supercooled liquid, this means we need 
to prevent the appearance of solid-like embryos above the 
critical size that would cause the liquid to freeze, which 
suggests we should choose the size of the largest solid- 
like embryo, Umax, as an order parameter to describe 
the state of the cluster f^. Furthermore, Umax seems to 
be a particularly appropriate order parameter in small 
nanoscale systems where the nucleation volume is suffi- 
ciently small that the appearance of a single post-critical 
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embyro leads to the termination of the metastable state 
throughout the entire system. When nmax = 0, the clus- 
ter is completely liquid but when n^ax = N, where N 
is the number of atoms in the cluster, then the cluster is 
completely solid, as a single crystal. The probability of 
finding the cluster in a given state is then 



Qij^max^ 



where 



Q{n^ax)^{l/NIA'^)J2' 



-E(nrr,a^}/kT 
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is the canonical partition function for the system con- 
strained to contain at least one largest embryo of size 
n-max so the sum is over all states characterised with a 
given rimax, E{nmax) is the energy, k is Boltzmann's con- 
stant, T is the temperature and A is the de Broglie wave- 
length. P{nmax) can be calculated by simulation and the 
free energy obtained from 



AF{n,-, 



-kT\nP(nr. 
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where /S.F{n,-aax) is the work required to take the entire 
system from a state where there is no solid-like cluster 
present, to a state where there is at least one largest clus- 
ter of size rimax- Eq- [3] closely resembles the intensive 
free energy introduced by ten Wolde and Frenkel 
to calculate the free energy barriers associated with nu- 
cleation, 

AF{n) = -kT\n (Pn/N) w -fcTln {N„,/N) , (4) 

where P„ is the probability of observing an n — sized 
embryo and -/V„ is the equilibrium number of embryos. 
AF{n) is the work of forming an n — sized embryo within 
the metastable phase. In the limit that embryos are rare 
(i.e. under conditions of mild undercooling) P(nmax) is 
approximately equal to the equilibrium number of em- 
bryos [m and the two free energies become equivalent 
within an additive constant, but it should be stressed 
that the two free energies are fundamentally different and 
that we would expect them to behave differently in deeply 
supercooled systems. 

Bagchi et al have recently used Eq.[3]to identify the 
liquid-gas spinodal in the supersaturated Lennard-Jones 
gas as the point at which AF{nmax) is a monotonically 
decreasing function of rimax- They find this occurs at 
a supersaturation consistent with previous estimates of 
the spinodal [isj and that the nucleation mechanism in 
the deeply metastable system changes from classical nu- 
cleation, characterised by fluctuating growth of a single 
large embryo, to a mechanism involving the coalescence 
of embryos. However, from the definition of P{nmax), it 
should be apparent that the free energy is an extensive 
quantity and it is likely that the location of the spinodal 



in a bulk system would shift depending on the number 
of particles in the simulation. In a small, finite-sized sys- 
tem, such as a liquid nanoparticle, the applicability of 
Eq. [3]is clearer. 

We have previously calculated AF{n) for a gold clus- 
ter with N — 456 atoms, for temperatures above T — 
690K but that work focused on the role of wet- 

ting phenomena and the location of the embryo at the 
nanoparticle surface. In the present paper, we calcu- 
late both AF{nmax) and AF{n) to lower temperatures 
for the same cluster in search of the limit of stability of 
the metastable liquid cluster using the same approach of 
combining umbrella Monte Carlo (MC) sampling simu- 
lation techniques with parallel tempering. We use the 
semi-empirical embedded-atom method (EAM) poten- 
tial [isl to describe the atomic interactions and study 
the cluster in the N, V, T ensemble with a simulation 
cell oi V ^ 1500A^ and periodic boundaries. At each 
temperature, we run eight parallel simulations or win- 
dows, each with a parabolic biasing potential w{nmax) — 
0.0005(n„ia2: — nof which biases the system to sample 
states where the largest embryo, Umax, in the cluster is 
around no- We choose no = 0, 10, 20, 30 . . . 70 and use 
T = 750, 730, 710, 690, 680, 670, 660, 650 for tempering. 

Our embryo criterion has been previously described in 
ref. [l3| and closely follows that developed by Frenkel [l3| 
to study crystal nucleation in hard sphere colloids. In 
brief, the criterion identifies which atoms in the cluster 
are solid-like by considering the degree to which the local 
order around two neighbouring atoms is correlated. If the 
local order of two atoms is highly correlated, then they 
are considered to be connected. If an atom is connected to 
half or more of its neighbours, then we consider the atom 
to be solid-like. Two solid-like atoms are considered to 
be in the same embryo if they are connected and Umax is 
the largest embryo. 

The embryo criterion is computationally expensive to 
apply so we use trajectories that consist of 10 normal 
MC moves per atom sampling the atomic interaction po- 
tential, followed by a test against w{nmax)- If the final 
move is rejected, the system is returned to state at the be- 
ginning of the trajectory. We attempt switches between 
neighboring no windows (T fixed) every 10 trajectories. 
We also attempt switches in neighboring temperatures 
(no fixed) every 10 trajectories, but these are offset with 
the Uo switches. These tempering switches have accep- 
tance ratios of about 0.4 and 0.6, respectively. The free 
energies in each window differ by an unknown additive 
constant, so the full free energy curve is constructed by 
fitting the curves to a polynomial in Umax fiof and a total 
of 1.74 X 10^ trajectories are sampled in each window. 

Fig. [1] shows the free energy calculated using Eq.[3l At 
temperatures just below the freezing temperature for the 
cluster, AF(n„iax) exhibits a minimum at small values of 
nmax before it increases to a maximum at larger embryo 
sizes. Umax = '^maa; dcnotes the critical embryo size. 



< -6- 



SltAAAAAAA'^'<,^«,«.««l«l««<«<^ 





750 K 




730 K 




710 K 




690 K 




680 K 




670 K 


• 


660 K 




650 K 



*•••. 




5 10 15 20 25 30 35 40 45 5 5 5 60 65 70 75 80 



FIG. 1: AF{nmax) as a function of Umax for temperatures in 
the range T = 750 - 650. 



Fluctuations in the cluster that keep the largest embryo 
below its critical size are locally stable and represent the 
configuration space available to the metastable liquid, 
while larger fluctuations cause the system to freeze. The 
critical size identifies the constraint required to keep the 
liquid in metastable equilibrium. 

As the temperature is lowered, n^^^^ decreases in size 
and the barrier becomes smaller. Eventually we reach a 
point, at T = 660K, where the barrier disappears and all 
fluctuations which increase the size of the largest clus- 
ter lower the free energy, suggesting we have reached the 
limit of stability for the fluid phase. Further decreases in 
T simply increase the thermodynamic driving force to- 
wards forming the solid as the free energy curve becomes 
steeper. 

Fig. [2] shows the two free energies calculated from 
Eqs. [3] and [3] where AF{n„iax) has been shifted verti- 
cally to maximise the overlap between the two curves. 
At T = 750K (see insert), the two free energies are iden- 
tical for embryo sizes larger than about 15 since there is 
generally just one large embryo in the system. The min- 
imum in AF{nmax) suggests the cluster usually contains 
a largest cluster of n « 5, but since AF{n) continues 
to decrease, there must be a larger number of smaller 
embryos present. At the spinodal temperature, the two 
curves are very different and only overlap at the largest 
embryo sizes. 

If we define the height of the barrier, AF(nJ„^^), as 
the difference in free energy between the maximum and 
the small embryo minimum, we can compare this with 
the usual nucleation barrier, AF{n*). Fig. [3^ shows that 
as the AF(n^£j2,) goes to zero at the spinodal, AF{n*) 
plateaus as a function of temperature at around lOkT. 
At the same time, the size of the critical embryo for both 
free energies decreases as a function of temperature. At 
the spinodal, AF{nmax) exhibits a flat region, where the 



FIG. 2: Comparision of AF{n) and AF{n,nax) at T 
and T = 750K (insert). 
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FIG. 3: (a) The height of the free energy barrier obtained 
from the free energy defined in Eq. [3] (squares) compared to 
that obtained from a free energy based on the equilibrium dis- 
tribution of embryos (circles) . See ref. [ij) for details, (b) The 
size of the critical embryo obtained using the two methods. 
Symbols are the same as (a). 



embryo sizes in the range Umax = 5 — 25 have approx- 
imately the same free energy, so we can expect consid- 
erable fluctuations in the embryo size. Nevertheless, the 
Timax remains finite (Fig.[3]3). This is in direct contrast to 
the predictions of mean-field theory 0,0], but our results 
are consistent with those of Pan et al '511 and Bagachi et 
al [ij. 

The rate at which clusters freeze can be determined 
by considering an ensemble of temperature quenched, 
molecular dynamics (MD) simulations . The liquid 
cluster is initially equilibrated at T = 900K, well above 
the freezing temperature, before the temperature is in- 
stantaneously quenched below freezing by rescaling the 
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particle velocities. The MD trajectory is then followed 
as the cluster freezes. Assuming this process is described 
by a first order rate law, the nucleation rate, J, can be 
obtained from the relation 



\n[R{t)] = -JVc{t - io) , 



(5) 



where R{t) is the fraction of un-nucleated clusters at time 
t, Vc is the volume of the cluster and to is the nucle- 
ation lag time, which is the time required to reach the 
steady state concentration of precritical embryos after 
the quench. To make use of Eq. \E\ we consider a cluster 
to have nucleated when rimax is greater than 85 for the 
last time during the simulation, which runs for 500 pi- 
coseconds. The nucleation size is defined as 85 because 
it is larger than the critical embryo size at all tempera- 
tures studied. A total of 300 quenched simulations are 
used at each temperature and even at the slowest rates 
(highest temperatures), less than 5% of the clusters re- 
mained un-nucleated by the end of the simulation. The 
volume of the cluster was determined using a "rolling 
sphere" algorithm [3] which defines the surface of a clus- 
ter using a hard sphere probe. In our case, the radius of 
the probe sphere and the gold atoms was taken to be 
1.5A. At T = 750K, Vc = 7 X 10^ ± 250A^, which is 12% 
smaller than would be predicted based on the volume per 
molecule of bulk Hquid EAM gold [13] ■ 

Fig. [5] shows that the nucleation rate increases as the 
cluster is quenched to lower temperatures. For tempera- 
tures below TOOK, our rates are approximately the same 
as those obtained by Bartell et al ^3], who used the 
same technique, but a larger cluster volume and a dif- 
ferent nucleation criterion. Around T = TOOK we see 
an unexpected increase in the rate with the slope dJ/ dT 
becoming more negative. Classical nucleation theory ex- 
presses the rate of nucleation as 



J = A'cxp[-AF7fcr] , 



(6) 



where the kinetic prefactor is given by AT = 
24pnZDn*^/'^/\, D is the diffusion coefficient, p„ is the 
number density of particles, A is the typical distance a 
particle must diffuse in order to join the embryo and 
Z = (|A^|/67rfcT7i*)i/2 is the Zeldovich factor. is 
the difference in chemical potential between the nucleat- 
ing stable and metastable phases. The temperature de- 
pendent parameters in the rate should vary continuously 
as a function of temperature and cannot account for the 
rapid increase in rate while Fig. [3^ suggests that the tem- 
perature dependence of AF{n*)/kT would cause the rate 
to slow, rather than accelerate. However, at T = TOOK, 
the barrier defined by AF{n^^^)/kT is in the order of 
kT, which suggests the observed deviation in the temper- 
ature dependence of the rate might be associated with a 
crossover from a barrier dominated nucleation process to 
a barrierless one. Consequently, both our direct barrier 
calculations and the independent MD rate calculations 
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FIG. 4: The nucleation rate as a function of temperature. 
Insert. The same rate data on a log scale. 



point to the strong possibility of a spinodal signifying 
the limit of stability of the fluid phase. 
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